Multifractal Analysis of Human Retinal Vessels 
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In this work it is shown that vascular structures of the human retina represent geometrical mul- 
tifractals, characterized by a hierarchy of exponents rather then a single fractal dimension. A 
number of retinal images from the STARE database (www.parl.clemson.edu/stare) are analyzed, 
corresponding to both normal and pathological states of the retina. In all studied cases a clearly 
multifractal behavior is observed, where capacity dimension is always found to be smaller then the 
information dimension, which is in turn always smaller then the correlation dimension, all the three 
being significantly lower then the DLA (Diffusion Limited Aggregation) fractal dimension. We also 
observe a tendency of images corresponding to the pathological states of the retina to have lower 
generalized dimensions and a shifted spectrum range, in comparison with the normal cases. 

PACS numbers: 05.40.-a, 61.43.Hv, 87.57.-s, 87.57.Nk 
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Over the past decade, there have been several attempts 
000000 in the direction of employing the frac- 
tal dimension as a measure for quantifying the "state" of 
human retinal vessel structures (considered as geometri- 
cal objects), with the expectation that such analysis may 
contribute to automatic detection of pathological cases, 
and therefore to computerization of the diagnostic pro- 
cess. While this is certainly a valid question with pos- 
sibly high impact on real world diagnostic issues, there 
are some issues that should be addressed before such in- 
vestigations may prove useful for the standard clinical 
practice. First, the fact that retinal vessels represent "fi- 
nite size" realizations of a fractal growth process, imposes 
questions about how exactly should one go about mea- 
suring the fractal dimension of a particular instance (e.g. 
an electronic image of a retinal vessel structure, taken 
from a given angle, with a given resolution and light- 
ning conditions). A related question is to what extent 
these calculations may correspond to the limiting fractal 
(which would have been attained if the growth process 
could have been extended to infinity), or equivalently, to 
what degree they may be compared with calculations on 
other such finite instances. Although various issues re- 
lated to these questions have already been addressed (for 
a current review see e.g. @), it seems that many of them 
remain open for further investigation. Second, some of 
these works 0, address the point that the retinal ves- 
sels may have different properties in different regions, and 
do indeed find different characteristics depending on the 
locale of measurement, although the procedures adopted 
in these works are only remotely related to established 
concepts of multifractality, and the corresponding com- 
monly acce pted pr ocedures for its measurement (see e.g. 
0M Holllll Il2| and references therein) . 

In this work we concentrate on the latter of the above 
issues, that is, we show that the human retinal vessel 



structures are geometrical multifractals, characterized by 
a hierarchy of exponents rather then a single fractal di- 
mension. We analyze a number of retinal images from 
the STARE database corresponding to both nor- 
mal and pathological states of the retina. In all cases 
we find clearly multifractal behavior. The capacity (or 
box counting) dimension is always found to be smaller 
then the information (or Shannon) dimension, which is 
in turn always smaller then the correlation dimension. In 
all cases the observed values of the capacity dimension 
were significantly lower then the DLA (Diffusion Limited 
Aggregation) fractal dimension, which has been consid- 
ered in earlier works 0, 0] as the primary model rele- 
vant for the phenomenon at hand. It is also found that 
images corresponding to pathological cases tend to have 
lower generalized dimensions, as well as a shifted spec- 
trum range, in comparison with the normal cases. 

In contrast to regular fractals (or monofractals), mul- 
tifractals are characterized by a hierarchy of exponents, 
rather then a single fractal dimension. A well known ex- 
ample of multifractality is the growth probability distri- 
bution during the DLA growth process, which has been 
shown to exhibit multifractal scaling |l3. ITR ITfiL Il7j . 
Geometrical (or mass) multifractals represent a special 
case when the measure of interest is homogeneously dis- 
tributed over the observed structure, so that only the 
number of particles (Lebesgue measure) contributes to 
the measure within a given region of the fractal 0, 0| • 
Considering a structure with mass (number of pixels) 
Mq and linear size L, covered with a grid of boxes of 
linear size £, the generalized dimension D q for the mass 
distribution is defined by 
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where Mi is the mass (number of pixels) within the i- 
th box, and q is a continuous (adjustable) variable that 
makes it possible to single out fractal properties of the ob- 
ject at different scales (equivalent of inverse temperature 
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in thermodynamics). The generalized dimensions Dq, Di 
and D 2 correspond to the capacity (or box-counting) di- 
mension, information (or Shannon) dimension, and cor- 
relation dimension, respectively. Finally, -D_oo and Doo 
represent the limits of the generalized dimension spec- 
trum (for monofractals, all the generalized dimensions 
coincide, being equal to the unique fractal dimension). 

It turns out that the direct application of in prac- 
tice is hindered by the fact that for q < the boxes 
that contain a small number of particles (because they 
barely overlap with the cluster) give anomalously large 
contribution to the sum on the left hand side of Q. To 
alleviate this problem and perform the multifractal anal- 
ysis of the retinal vessel structures, we adopt the gen- 
eralized sand box method 0, Ej, which has been suc- 
cessfully used do demonstrate geometric multifractality 
of DL A 9] . This procedure consists in randomly select- 
ing iV points belonging to the structure, and counting for 
each such point i the number of pixels Mi(R) that belong 
to the structure, inside boxes of growing linear dimension 
R, centered at the selected pixels. The left hand side of 
equation can now be interpreted as the average of the 
quantity (Mi(R)/M ) q ~ according to probability distri- 
bution Mi(R) /Mq. When the box centers are chosen ran- 
domly, the averaging should be made over the chosen set, 
and the equivalent of becomes 

((^rHr" < 2 » 

The practical advantage of this method is that the boxes 
are centered on the structure, so that by construction 
there are no boxes with too few particles (pixels) inside. 

To verify whether human retinal vessel structures 
demonstrate geometrical multifractal scaling properties, 
we have used a set of forty retinal images from the 
STARE database 0], manually segmented by two dif- 
ferent observers (herefrom referred to by initials AH and 
VK as in 13]) from twenty original retinal scans (con- 
taining ten normal and ten pathological cases), for the 
purpose of studies on automatic image segmentation and 
diagnostics • The images segmented by observers AH 
and VK differ in level of detail, and the resulting set, 
totaling forty segmented images, is available for down- 
load from the STARE project [13| in ppm file format, 
with resolution of 700x605 pixels. As recently it has been 
argued |19| that fractal analysis may be more sensitive 
to changes in vascular patterns when skeletal images of 
vascular trees are considered, rather then the original 
segmented images (which contain the vessel width in- 
formation), in order to verify whether the vessel width 
information indeed does exert influence on the multifrac- 
tal analysis, we have also performed skeletonization of the 
two downloaded sets using the eight connected Rosenfeld 
algorithm J2j}, to produce two additional sets of twenty 
images each. A typical normal and a pathological im- 
age, segmented by observers AH and VK (where images 
segmented by observer VK demonstrate a substantially 




FIG. 1: Image of a typical normal retinal vessel structure 
(image files im0162.ah.ppm from the STARE database 
segmented by a) observer AH and b) observer VK, together 
with their skeletonized versions c) and d) , respectively, and a 
typical pathological structure (image file im0001.ah.ppm, di- 
agnosed with Background Diabetic Retinopathy), segmented 
by observer AH. 

higher level of detail), respectively, together with their 
skeletonized versions using the Rosenfeld algorithm, are 
shown in Fig. 1. 

For all of the four sets (containing twenty images 
each), we have performed measurements (calculations) 
according to J2J, selecting 1000 random points on each 
structure (typical structure size Mq is approximately 
30000 pixels, and the typical linear size L is 600 pix- 
els) , and counting number Mi of pixels inside boxes cen- 
tered at selected points. These numbers were then used 
to extract generalized dimension D q , for different val- 
ues of q (—10 < q < 10), as slopes of the lines ob- 
tained through regression (minimum squares fitting) of 

log([M(R)/M ] q ~ 1 \/(q- 1), as a function of log(R/L). 

The whole procedure was repeated 100 times (with dif- 
ferent random choices of the 1000 random points), for 
each image, and for each value of q. The final values of 
D q were calculated as averages over these repetitions. 

A word is due on calculations for the special case 
q = 1, corresponding to information dimension D±. As 
the above formulas are non-analytic for q = 1, we per- 
form calculations at q — lie, for e = 0.001, and assuming 
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FIG. 2: The generalized dimension spectrum, D q versus q, 
for a typical normal retinal image (image file imOf 62.ah.ppm 
fT^n . averaged over 100 random choices of 1000 points each 
(see text for details). The error bars indicate the largest 
and smallest values encountered within the 100 runs, and the 
curve connecting the points serves as a guide to the eye. The 
points corresponding to the capacity dimension Do — 1.647, 
the information dimension Di = 1.594 and the correlation di- 
mension D2 = 1.552 are represented by full circles, while the 
dotted lines serve to emphasize their position. 



linearity of the function D{q) in this (short) interval, we 
interpolate D\ » (-Di- e + fi+ e ) /2 (in fact, the differ- 
ence between the values of D q calculated on both sides 
of q = 1 was found to be only slightly larger then the 
statistical fluctuations induced by random choice of the 
set of measurement points on the structure). 

Results of a typical calculation are shown in Fig. 2, 
where it is seen that the observed retinal vessel structure 
clearly demonstrates multifractal scaling, rather then be- 
ing a simple monofractal (there is a significant differ- 
ence between generalized dimensions). In particular, the 
capacity dimension Dq, the information dimension Di 
and the correlation dimension Di are all different, sat- 
isfying Do > D\ > Z?2- Also, all the three values re- 
main substantially lower then the DLA fractal dimen- 
sion estimate, commonly accepted in the literature, of 
D q= 2 — 1.71 (which is in fact underestimated by com- 
monly used methods) in contrast with previous find- 
ings Uli. 

Numerical corresonding to Fig. 2 (for the set of twenty 
images from the STARE database segmented by ob- 
server AH) are given in Tab. [I] The first column lists 
the image names, while the second column indicates im- 
age classification status as "Pathological" or "Normal" 
|2l| . The values of generalized dimensions D q are given 
for q = —10, 0,1, 2, 10, where as already mentioned Do, 
D\ and Di correspond to the capacity, information and 
correlation dimension, respectively, while -D-10 and D±q 
indicate the range of the general dimension spectrum. 
It is seen from Tab. H] that all of the values calculated 
for the capacity dimension (which corresponds to box 



TABLE I: Generalized dimensions D q for < 
for the twenty analyzed images from the 
The second column indicates classification 
the images (pathological and normal). 

Image Status D_io Do Di 



j = -10, 0,1, 2, 10, 
STARE database, 
status for each of 

D 2 D 10 



imOOOl.ah 
im0002.ah 
im0003.ah 
im0004.ah 
im0005.ah 
im0044.ah 
im0077.ah 
im0081.ah 
im0082.ah 
im0139.ah 
im0162.ah 
im0163.ah 
im0235.ah 
im0236.ah 
im0239.ah 
im0240.ah 
im0255.ah 
im0291.ah 
im0319.ah 
im0324.ah 



P 
P 
P 
P 
P 
P 
N 
N 
N 
P 
N 
N 
N 
N 
N 
N 
N 
P 
P 
P 



1.968 
1.930 
1.877 
1.777 
1.900 
1.886 
1.911 
1.917 
1.981 
1.904 
1.968 
1.998 
1.945 
1.868 
1.945 
1.918 
1.944 
1.819 
1.703 
1.923 



540 
548 
509 
522 
589 
541 
576 
555 
578 
565 
1.647 
1.642 
1.597 
1.584 
1.587 
1.593 
1.633 
1.516 
1.443 
1.567 



1.494 
1.498 
1.469 
1.492 
1.560 
1.493 
1.528 
1.514 
1.518 



.516 
.594 
.587 
.548 
.544 
.549 
.564 
1.604 
1.482 
1.409 
1.520 



1.462 
1.460 
1.443 
1.465 
1.538 
1.459 
1.496 
1.487 
1.476 
1.481 



.552 
.550 
.514 
.514 
.520 
.543 
.583 
1.454 
1.382 
1.486 



1.361 
1.370 
1.380 
1.367 
1.474 
1.363 
1.426 
1.421 
1.404 
1.413 
1.459 
1.476 
1.442 
1.448 
1.437 
1.494 
1.521 
1.348 
1.299 
1.399 



counting method), and indeed the correlation dimension 
(corresponding to methods such as radius of gyration or 
the density-density correlation function) , are significantly 
lower then the DLA fractal dimension D q= i ~ 1.71 
Therefore, our results show that retinal vessel structures 
are geometrical multifractals, and that the overall fractal 
dimension is lower then that of the DLA. 

Results of the multifractal analysis for the other three 
sets of images (STARE database images segmented by 
observer VK, and the skeletonized versions of AH and 
VK) all yield qualitatively similar results, all of them 
clearly demonstrating multifractal behavior. In Table ILTI 
we present the results for the capacity (box counting) 
dimension Dq, for all of the four sets of images. It fol- 
lows from Table [Q] that the process of skeletonization 
(removal of vessel width information from the image) 
slightly reduces the fractal dimension, while this effect 
is much weaker in comparison with the effect of the level 
of detail present in the segmentation process, as found 
between the two current observers. However, when the 
results are compared consistently within each group sep- 
arately, the mean fractal dimension is found to be lower 
for the pathological images then for the normal cases, for 
all of the four groups. Although this finding can hardly 
be considered conclusive from the statistical viewpoint, 
it is nevertheless encouraging from the point of view that 
fractal spectrum analysis could be employed for quantifi- 
cation of the retinal vessel state, in order to contribute 
to automatic diagnostics. To this end, far more detailed 
studies of images corresponding to specific deseases and 
normal cases, are needed. Assuming that each of the ob- 
servers consistently applied his own criteria in segmenta- 
tion, it follows that the fractal dimension results may be 
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TABLE II: Capacity(or box counting) dimension Do for the 
two sets of images from the STARE database segmented by 
observers AH and VK, together with their skeletonzied ver- 
sions. The second column indicates classification status for 
each of the images (pathological and normal), and the last 
three lines present averages for the pathological, normal and 
all images, respectively. 

Image Status AH AH-S VK VK-S 



imOOOl 


P 


1.540 


1.545 


1 


583 


1 


593 


im0002 


P 


1.548 


1.524 


1 


574 


1 


568 


im0003 


P 


1.509 


1.500 


1 


593 


1 


608 


im0004 


P 


1.522 


1.508 


1 


573 


1 


598 


im0005 


P 


1.589 


1.554 


1 


680 


1 


663 


im0044 


P 


1.541 


1.538 


1 


668 


1 


661 


im0077 


N 


1.576 


1.591 


1 


658 


1 


662 


im0081 


N 


1.555 


1.551 


1 


668 


1 


671 


im0082 


N 


1.578 


1.585 


1 


665 


1 


680 


im0139 


P 


1.565 


1.564 


1 


679 


1 


678 


im0162 


N 


1.647 


1.638 


1 


714 


1 


700 


im0163 


N 


1.642 


1.612 


1 


684 


1 


646 


im0235 


N 


1.597 


1.588 


1 


685 


1 


675 


im0236 


N 


1.584 


1.581 


1 


658 


1 


662 


im0239 


N 


1.587 


1.597 


1 


655 


1 


655 


im0240 


N 


1.593 


1.563 


1 


677 


1 


663 


im0255 


N 


1.633 


1.634 


1 


696 


1 


693 


im0291 


P 


1.516 


1.491 


1 


604 


1 


578 


im0319 


P 


1.443 


1.446 


1 


555 


1 


561 


im0324 


P 


1.567 


1.503 


1 


642 


1 


617 


Average 


P 


1.534 


1.517 


1.615 


1.612 




N 


1.599 


1.594 


1.676 


1.671 




All 


1.567 


1.556 


1 


646 


1.642 



compared only between images segmented by the same 
observer, either skeletonized or not, but should be nor- 
malized before making comparisons of results from dif- 
ferent groups. 

When addressing multifractality, numerous works deal 
with the so-called f(a) spectrum (see e.g. 0,0,112 an d 
references therein), where 



N(a) = L- fia \ 



(3) 



represents the number of boxes N(a) where the proba- 
bility Pi of finding a particle (pixel) within a given region 
i scales as 



Pi = L ai 



(4) 



and /(a) may be understood as the fractal dimension 
of the union of regions with singularity strenghts be- 
tween a and a + da. The exponent a takes values 
from the interval [—00,00], and the function /(a) is 
usually a single humped function with a maximum at 
df (a(q)) I da(q) = 0. The relationship between the D(q) 
spectrum and the /(cv) spectrum is made via the Legen- 
dre transform 



where 



/(a (<?)) =qa (q) - r (q) 
dr(q) 



a{q) 



dq 



(5) 
(0) 



f(a) 




FIG. 3: The f(a ) sp ectrum for the twenty images from the 
STARE database jl3J, segmented by observer AH. Curves cor- 
responding to normal retinal images are represented by open 
circles, and those corresponding to pathological images |2l| 
are represented by full symbols. It is seen that pathological 
image curves tend to be shifted to the lower a range and have 
lower maxima, in comparison with the normal images (see 
text for more details). 



and 



r{q) = {q- l)D q 



(7) 



is the mass correlation exponent of the q th order. To 
calculate the derivatives in @, we have performed cal- 
culations at pairs of points q and q + e with e = 0.001, 
so that derivatives were calculated as dr(q)/dq « (r(q + 
e) — r(q))/e, except at point q = 1, where we have used 
dr(q)/dq^(r(l + e)-T(l-e))/(2e). 

In Fig. |3| we show detailed results of our calculations, 
performed on the STARE database images segmented by 
observer AH, with respect to the f(a) spectrum. While 
the current set of images is not particularly adequate for 
testing the effects of a given type of pathology (there 
are only ten normal images, and ten pathological images 
affected by not necessarily the same disease, see [2l|), it 
is seen that pathological case images tend to have lower 
maxima, occasionally more narrow spectrum range, and 
a shift in the spectrum position, in comparison with the 
normal cases. 

Finally, in Fig. 0] we present results of the f(a) spec- 
trum averaged separately for the normal and the patho- 
logical images for all of the four sets, where it is seen 
that the previous observation holds for both observers, 
independent of skeletonization. The skeletonized images 
present more narrow f(a) spectrum then the original seg- 
mented images (which contain the vessel width informa- 
tion) for both observers, which may explain the conclu- 
sion of that fractal analysis after skeletonization may 
be more sensitive to changes in vascular patterns. More 
precisely, since monofractals have infinitely narrow /(a) 
spectrum (a single fractal dimension), the above results 



5 




FIG. 4: The f(a ) sp ectrum for the twenty images from the 
STARE database (l3j ]. segmented by observer AH. Curves cor- 
responding to normal retinal images are represented by open 
circles, and those corresponding to pathological images |2l| 
are represented by full symbols. It is seen that pathological 
image curves tend to be shifted to the lower a range and have 
lower maxima, in comparison with the normal images (see 
text for more details). 

show that skeletonized structures may be more closely 
approximated as monofractals (when a single dimension 
is calculated rather then the whole spectrum). As the 
general properties of the spectrum are preserved through 
skeletonizatin, another advantage of using such images 
may be considered the fact that they contain far fewer 
pixels, and therefore the calculations require less com- 



puter time. 

The results of calculations of the f(a) spectrum pre- 
sented in Figs. 13141 again may be considered encouraging 
from the point of view of the objective of turning the 
diagnostic process automatic, although further more de- 
tailed studies are necessary to determine their statisti- 
cal significance, and whether the observed differences in 
multifractal scaling behavior may be exploited for dis- 
cerning normal images from images with certain types of 
pathologies. More precisely, the current work is primar- 
ily concerned with establishing the fact that retinal vessel 
images represent geometrical multifractals, nevertheless, 
our calculations suggest that there may be grounds for 
automatic differentiating between normal images and cer- 
tain pathological cases. 

In conclusion, we show in this work that vascular struc- 
tures of the human retina represent geometrical multi- 
fractals, characterized by a hierarchy of exponents, rather 
then a single fractal dimension. We analyze twenty reti- 
nal images from the STARE database [R| , where half of 
the images correspond to normal states of the retina, and 
half to different pathological states [2l|, together with 
their skeletonized versions. In all studied cases we find 
clearly multifractal behavior, with capacity dimension 
considerably lower then the DLA value. We also observe 
a tendency of normal images of having higher general- 
ized dimensions and a shift of the f(a) spectrum range 
towards higher singularity strength values, in comparison 
with the pathological cases. While the last observations 
are hardly conclusive from a statistical standpoint, they 
may prove relevant in the quest of automatic diagnostic 
procedures. 
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